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Thermal structure and cooling of neutron stars 
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Abstract. The thermal structure of neutron stars with magnetized envelopes is studied using modern physics 
input. The relation between the internal (Tint) and local surface temperatures is calculated and fitted by analytic 
expressions for magnetic field strengths B from to 10 16 G and arbitrary inclination of the field lines to the 
surface. The luminosity of a neutron star with dipole magnetic field is calculated and fitted as a function of B, 
Tint, stellar mass and radius. In addition, we simulate cooling of neutron stars with magnetized envelopes. In 
particular, we analyse ultramagnetized envelopes of magnetars and also the effects of the magnetic field of the 
Vela pulsar on the determination of critical temperatures of neutron and proton superfluids in its core. 
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1. Introduction 

It is well known that theoretical models of cooling of iso- 
lated neutron stars (NSs) depend on the poorly known 
equation of state (EOS) of superdense matter in the NS in- 
teriors. Comparing calculated cooling curves (decrease of 
the NS effective surface temperature T e in time) with ob- 
servations gives a potentially powerful method of testing 
microscopic theories of superdense matter (e.g., Pethick 
|1992; Page 1997, 1998). Its practical implementation is 



restricted by the accuracy of observational determination 
of T c and NS ages and by the quality of the physics in- 
put used in calculations. Great observational progress has 
been achieved recently after the launch of Chandra and 
Newton X-ray observatories. In the present paper we up- 
date theoretical models of thermal structure and evolution 
of NSs with magnetized envelopes. 

Most of the observed NSs p ossess magnetic fields B ~ 
10 12 -10 13 G (e.g., Taylor et al. |199$) . 
sibly m agneta rs, with B > 10 14 G 

Mereghetti 



Duncan 1995: Kouvcliotou ct al. 1998 



Some NSs are pos- 
e .g., T hompson & 
" 199S| ; 



2001 ). The internal NS magnetic field can be even higher. 
The strong magnetic held affects physical properties of all 
NS layers in many ways. For instance, the field B ~ 10 16 
G may affect the neutrino emissivity in the NS core (e.g., 
Baiko & Yakovlev[T999|). The field B k, 10 12 G may change 



the electrical resistivity of NS cores, accelerating evolution 
of the internal held. In principle the thermal evolution may 



be co upled to the magnetic one (e.g., Urpin & Shalybkov 
1995b . 

Thus modelling of the thermal structure and evolu- 
tion of the magnetized NSs is a complicated task. In this 
paper we focus on two problems. First we consider the 
thermal structure of the outer magnetized NS envelope of 
density p < pb = 4 x 10 11 g cm" 3 with the magnetic held 
B 10 16 G. These envelopes produce thermal insulation 
(blanketing) of NS interiors. We solve this problem using 
updated physics input described in Sect. ^[ The solution 
(Sect. |!|) relates the internal NS temperature T[ nt with the 
local effective surface temperature T s and, therefore, with 
the integrated NS luminosity (or the mean effective sur- 
face temperature, T c ). Second, in order to illustrate the 
obtained Ts(Tj nt ) relation, we simulate (Sect. |j) cooling of 
NSs which possess a given dipole or radial magnetic field 
in the outer envelopes. 

The effects of a strong magnetic field on thermody- 
namic and kinetic properties of the outer NS layers have 
been reviewed, for instance, by Yakovlev & Kaminker 
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(|199J) and Ventura & Potekhin ( gOOjJ ). The field affects 
the properties of all plasma components, the electron com- 
ponent usually being affected most strongly. Motion of 
free electrons perpendicular to the field lines is quantized 
in Landau orbitals with a characteristic transverse scale 
equal to the magnetic length a m = (Hc/eB) 1 ^ 2 = A c /V&, 
where A c = h/(m e c) is the electron Compton wavelength, 
b = hoj c /m e c 2 = i?i2/44.14 is the magnetic field strength 
expressed in relativistic units, uj c = eB/m c c is the electron 
cyclotron frequency, and Bi 2 = B/(10 12 G). Except for 
the outermost parts of the NS envelopes, the electrons con- 
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stitute degenerate, almost ideal gas. The electron entropy, 
magnetization, thermal and electrical conductivities, and 
other quantities exhibit quantum oscillations of the de 
Haas-van Alphen type. The oscillations occur under vari- 
ations of density p or B whenever the electron Fermi mo- 
mentum reaches the characteristic values \j2nb m c c asso- 
ciated with occupation of Landau levels n — 1,2,... These 
oscillations appreciably change the properties of degener- 
ate electrons in the limit of a strongly quantizing field (e.g., 
Yakovlev & Kaminker 1994) in which almost all electrons 



populate the ground Landau level. The latter case takes 
place at temperature T <C Tb an d density p < ps, where 



p B = m u n B A/Z ^ 7045 B% 2 (A/Z) gem 



,-3 



T B = hLJ s /k B ~ 1-343 x 10 8 {B 12 /y/l + x 2 ) K. 



(1) 
(2) 



Here, m u = 1.66 54 x 10 24 g is the atomic mass unit, A 
and Z are the mean ion mass and charge numbers, uj g = 
cj c / -Jl + x% is the electron gyrofrequency, 



^V) 173 . LQQ9 (, 6W /3 



(3) 



is the non-magnetic relativity parameter, p^ = 
p/(10 6 gcm~ 3 ), and n B = l/fir^ai) is the elec- 
tron number density at which the Fermi energy reaches 
the first excited Landau level. In the case of a weakly 
quantizing magnetic field (p ^ pb, T Tb), the 
oscillations of electron quantities occur around the 
classical non-magnetic values and are typically not very 
pronounced. In the case of T ^> Tb, the field can be 
treated as classical (non- quantizing) and the oscillations 
disappear. 

The thermal structure of magnetized NS envelopes 
has been studied by a number of authors mainly us- 
ing a plane-parallel approximation. More attention has 
been paid to the case of the radial magnetic field (nor- 
mal to the surface). It has been thoroughly considered 
by Hernquist ( |l985j), Van Riper ( [1988D , Schaaf ( |l990a| ), 
Heyl & Hernquist ( |l998b|) for B < 10 14 G. In another pa- 
per Heyl & Hernquist (1998a) analysed the case of higher 
fields, B ~ 10 15 -10 16 G. One can consult the cited papers 
for references to earlier works. The principal conclusion 
of these studies is that the magnetic field reduces ther- 
mal insulation of the blanketing envelope by increasing 
the longitudinal (along the field lines) thermal conductiv- 
ity of degenerate electrons due to Landau quantization of 
electron motion. 

The thermal structure of the envelope with the mag- 
netic field tangential to the surface has been analysed by 
Hernqu ist ( 1985 ), Schaaf ( 1990a ), and Heyl & Hernquist 
( |l998a| ) for B ^ 10 14 G. In this case the field increases 
thermal insulation of the blanketing envelope due to the 
classical effect of reduction of the electron thermal con- 
ductivity transverse to the field because of the Larmor 
rotation. 

The case of arbitrary inclination of the field to the 
surface was studied by Greenstein & Hartke (1983) in the 



approximation of constant (density and temperature inde- 
pendent) longitudinal and transverse thermal conductiv- 
ities. The authors proposed a very simple formula (Sect. 
[j.3| ) which relates the local surface and internal stellar 
temperatures, T s and Xi nt . It is constructed from two 
7s(Ti n t) relations obtained for the radial and tangential 
magnetic fields. Page (1995) presented arguments that the 
formula of Greenstein & Hartke is valid also for realis- 
tic, variable thermal conductivities. If so, one immediately 
gets the required T s (Ti nt ) relation for any magnetic field 
inclination using the realistic relations for the radial and 
tangential fields. This method has been used by several 
authors (e.g., Shibanov & Yakovlev 1996| ). Recently, the 
case of arbitrary field inclination has been studied also 
by Heyl & Hernquist ( 1998b ). In Sect. || we reconsider 
the thermal structure of the blanketing envelopes for any 
magnetic field inclination and compare our results with 
those of earlier studies. 

Early simulations of cooling of magnetized NSs were 
performed assuming the radial magnetic field everywhere 
over t he ste llar surface (e.g., Nomoto & Tsuruta 1987; Van 
Riper |l99l|) . Since the radial magnetic field reduces the 
thermal insulation, these theories predicted acceleration 
of cooling of the magnetized NS accompanied by enhanced 
NS luminosity at the early (neutrino-dominate d) co oling 
stage (at stellar ages t ^ 1 4 — 10 5 yr). Page ( 1995 ) and 
Shibanov & Yakovlev ( 1996| ) simulated cooling of NSs with 
dipole magnetic fields. They showed that the decrease of 
thermal insulation of the stellar envelope near the mag- 
netic equator was partly compensated by its increase near 
the pole, and the magnetic field did not necessarily accel- 
erate the cooling, in agreement with an earlier conjecture 
of Hernquist fll985|). In a series of papers Heyl & Hernquist 
( |!997a| , |l997b| |l998a| , |l998bQ proposed simplified models 
of cooling of magnetized NSs including the cases of ultra- 
high surface magnetic fields, B ~ 10 15 — 10 16 G. 

We illustrate our new models of magnetized NS en- 
velopes with simulation of cooling of NSs (Sect. |I]). We 
briefly discuss cooling of ultramagnetized NSs as well as 
cooling models of the Vela pulsar with the dipole magnetic 
field and superfluid core. 



2. Physics input 

2.1. Equations of thermal structure and evolution 

The internal hydrostatic structure of a NS can be regarded 
as temperature-independent (e.g., Shapiro & Teukolsky 



1983). It is conventional (Gudmundsson et al. 1983) to 



separate the calculation of heat transport in the NS inte- 
rior (at radius r < i?b) and in the outer heat-blanketing 
envelope (i?b < r < R, where R is the stellar radius). 
The choice of the boundary i?b must meet several require- 
ments. The blanketing envelope should be thin (R — i?b <C 
R) and contain negligibly small mass; there should be no 
large sources of energy generation or sink present there; 
it should serve as a good thermal blanket of the inter- 
nal region; its thermal relaxation time should be shorter 
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than time-scales of temperature variation in the internal 
region. In the studies of non-magnetized NSs, the bound- 
ary is usually taken at the density pb = 10 10 g cm" 3 (e., 
Gudmundsson et al. 



1983 



1997 - here- F = oF 



Potekhin et al 
after PCY). In the presence of a strong magnetic field we 
assume additionally that the temperature does not vary 
over the boundary r = i?b at any given moment. To make 
this requirement more realistic, we shift the boundary to 
the neutron drip density, pb = 4 x 10 11 gcm~ 3 . 



2.1.1. Internal region 

In the internal region (r < i?b), magnetic fields B ^ 
10 16 G are either non-quantizing or weakly quantizing 
(degenerate electrons populate many Landau levels). For 
simplicity, we neglect the effects of magnetic fields in this 
region and use a spherically symmetric temperature dis- 
tribution which obeys the classical equations of thermal 



evolution (Thorne 1977): 



-2<I> 



Airr 2 



2Gm d 



c 2 r dr 



(e»L) 



dT 

at' 



L 



4 7rr 2 



2Gm 



(4) 



(5) 



where Q y {r) is the neutrino emissivity, C v {r) is the heat 
capacity, n(r) is the thermal conductivity, L(r) is the lo- 
cal luminosity defined as the non-neutrino heat flux trans- 
ported through a sphere of radius r. Furthermore, G is the 
gravitational constant, m(r) is the gravitational mass in- 
side the sphere of radius r, and $(r) is the metric function 
determined by the stellar model. Notice that at the stellar 
surface = (1 - rg/R)- 1 / 2 , where r g = 2GM/c 2 = 

2.95(M/M Q ) km is the Schwarzschild radius defined by 
the total gravitational NS mass m(R) = M. On the right- 
hand side of Eq. ([|) we neglect the rate of heat produc- 
tion (e.g., by internal friction due to differential rotation, 
see Page 1997 , |l998| and references therein) . This effect is 
important for relatively old and cold NSs, which we will 
not consider here. To solve Eqs. (Q) and (pj), w e use the 



computer code described by Gnedin et al. (2001, hereafter 
GYP). After thermal relaxation within the NS, at t 100 
yr the redshifted temperature T(t) = T(r, i)e*( r ) becomes 
constant throughout the internal region. 

2.1.2. Blanketing envelope 

The thermal structure of the blanketing envelope is stud- 
ied in the stationary, local plane-parallel approximation 
to relate the local effective surface temperature T s to the 
temperature Tj nt at the inner boundary of the envelope. 
Everywhere on the surface, except possibly in tiny regions 
where the field is almost tangential, one may assume that 
a scale of temperature variation over the surface is much 
larger than the thickness of the blanketing envelope. This 



leads to the one- dimensional approximation for the heat 
diffusion equation: 

16cr T 3 dT 16crT 3 



dT 
~cL7 



dr 



3A> 



(6) 



where F is the local heat flux density (constant through- 
out a given local part of the envelope), k is an effective 
thermal conductivity along the normal to the surface, a 
is the Stefan-Boltzmann constant, K is the mean opac- 
ity (Sect. 2.3), r = J^^Kp dz is the optical depth, and 
z = (R— r) e*( fl ) is the local proper depth. 

Integration of Eq. (||) gives a temperature profile T w 
T s (|t + 5) 1 / 4 , and the integration constant corresponds 
to the Eddington approximation (r = | at the radiative 
surface, where T = T s ). A more accurate boundary con- 
dition requires solution of the radiative transfer equation 



in the NS atmosphere (e.g., Shibanov et al. 1998) 



At z <C R the general relativistic equation of hydro- 
static equilibrium can be reduced to the Newtonian form: 
dT/dz = gp, where g = GMj{R 2 ^J\ - r g /R) is the sur- 



face gravity. Together with Eq. 
structure equation 

dlogT 3 PK T s 4 



) , it leads to the thermal 



(7) 



dlogT 16 g 

where T is the pressure. 

The thermal conductivity tensor of magnetized plasma 
is anisotropic. It is characterized by the conductivities par- 
allel (km) and transverse (n±) to the field, and by the off- 
diagonal (Hall) component. In the plane-parallel approxi- 
mation Eq. (|]) contains the effective thermal conductivity 

K = K\\ cos 2 9 + k± sin 2 9, (8) 

where 9 is the angle between the field lines and the normal 
to the surface. In case 9 = the heat is transported solely 
by the parallel conductivity, n — while in case 9 = 90° 
it is transported by the transverse conductivity, K = kj_. 
These special cases will be referred to as the parallel (||) 
and transverse (_L) conduction cases. Actually, B and 9 
vary slowly over the surface. The total stellar luminosity 
is 



TdE 



cr / Tf dS = 47rcri? 2 T c 4 , 



(9) 



where dS = R 2 dfl is the surface element determined by 
corresponding solid angle d£l, and T e is the (mean) effec- 
tive temperature to be distinguished from the local effec- 
tive temperature T s . Naturally, T = T s for a non-magnetic 
NS. 

Assuming the dipole field, we can use the general- 
relativistic solution (Ginzburg & Ozernoy 1964) 



B(x) = B p y cos 2 x + a2 sm X> tan0 = atanx, (10) 
(1 - x) ln(l - x) + x - 0.5a; 2 



[ln(l - x) + x + 0.5x 2 ]t/T^x' 



(11) 



where B p is the field strength at the magnetic pole, x 1S 
the polar angle, and x = r g /R. 
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The quantities T s , T c and L refer to a local refer- 
ence frame at the NS surface. The redshifted ("apparent") 
quan tities as detected by a distant observer are (Thorne 
1977b : T s °° = T s y/l - r g /R, T e °° = T c y/1 - r g /R, = 
(l-r g /R)L. 

Crucial for the thermal evolution is the relation be- 
tween T s and Tj nt = T(i?b). For non-magnetic blanketing 
envelopes composed of iron, this relation was studied by 
Gudmundsson et al. ( 1983 ), while the non-magnetic en- 
velopes of various chemical compositions were analysed 
by PCY. Ventura & Potekhin fl200l|) presented an ana- 
lytic description of the Ti nt (T s ) relation for the envelopes 
without magnetic fields and with strong magnetic fields; 
Heyl & Hernquist (1998a) performed a semianalytic inves- 
tigation of the case of strong magnetic field. The results of 
Ventura & Potekhin fl2001j) and Heyl & Hernquist ( |l998a[ ) 
are rather approximate because of a number of simplifying 
assumptions discussed by Ventura & Potekhin ( |200l| ). 

The validity of the one-dimensional approximation can 
be checked by direct two-dimensional simulation of the 
heat transport in the blanketing envelope. Such simula- 
tion has been attempted by Schaaf ( |l990b[ ) for a homo- 
geneously magnetized NS under many simplified assump- 
tions, so that a more realistic study is required. The heat 
conduction from hotter to cooler zones along the sur- 
face or possible meridional and convective motions can 
smooth the temperature variations over the NS surface. 
Nevertheless, the one-dimensional approximation seems to 
be sufficient for our cooling calculations presented below. 



2.2. Equation of state 

In the inner region of a NS, the effects of magnetic field 
are assumed to be weak and we use the same EOSs as 
in GYP. Specifically, in the core, we use the moderately 
stiff phenomenological EOS of matter composed of neu- 
tron s, pro tons and electrons, as proposed by Prakash et 



al. ( 1988 ), in its simplified version suggested by Page & 
Applegate ( 1992 ). In the inner envelope, we use th e model 
of ground-state matter (Negele & Vautherin 197S ) and de- 
scribe the properties of atomic nuclei by the smooth com- 
position model (Kaminker et al. 1999 ). The core-crust in- 



terface is placed at p = 1.5 X 10 g cm 

The outer envelope is assumed to be composed of 
iron, which can be partially ionized at p ^ 10 6 gcm~ 3 . 
Following PCY, we employ the mean-ion approxima- 
tion and adjust an effective ion charge to a more elab- 
orate EOS: the Opacity Library (OPAL) EOS (Rogers 
et al. |1996[ ) at B = and the Thomas-Fermi EOS of 
Thorolfsson et al. (p^98|) at B = (10 10 -10 13 ) G. The tab- 



ular entries of these EOSs in p and T are interpolated in 
the same way as in PCY. When necessary, we interpolate 
the effective charge at intermediate B. Since no reliable 
EOS of iron has been published for B > 10 13 G, we use 
the effective charge obtained at B — 10 13 G for higher B. 

Thus, in our model of the outer envelope, the pres- 
sure is produced by magnetized Fermi gas of electrons 



and by the gas of classical ions with an effective charge. 
We neglect the anomalous magnetic moments of the nu- 
clei, which is a good approximation at B <C 10 16 G (cf. 
Broderick et al. |2000| ; Suh & Mathews p0f5l| ) . The electron 
pressure can be expressed through the chemical potential 
p and T (e.g., Blandford & Hernquist [1982p: 



P — P 



where 



br, 



3/2 n„ 



V2tt 2 



?i=0 



(2 - M(l + 2bn) 1 ^I 1/2 ( X n,r n ), (12) 



Vl + 2bn 



1 



T 



P T = m c c 2 /A 3 



1.4218 x 10 

~\9 



25 



\/l 

dyn cm 



2bn T T 
2 and 



T r = 



m e c /fes ~ 5.930 x 10 y K are the relativistic units of pres- 
sure and temperature, respectively. The standard Fermi- 
Dirac integra l is given by a fit presented in Chabrier 
& Potekhin ( 1998| ). The chemical potential p is found at 
given p, T, and B using an analytic fit and iterative pro- 
cedure described by Potekhin & Yakovlev ( |1996| ). 

2.3. Opacities 

The heat is carried through the NS envelope mainly by 
electrons at relatively high densities and by photons near 
the surface (e.g., Gudmundsson et al. 1983). In general, 



the two mechanisms work in parallel, hence 

K = Kr + K c , K- 1 =K- 1 +K- 1 , (13) 

where and K Y , K c denote the radiative (r) and elec- 

tron (c) components of the conductivity and opacity. 

Typically, the radiative conduction dominates [k t > 
k c ) in the outermost non-degenerate layers of a NS, 
whereas the electron conduction dominates (k c > k t ) in 
deeper, moderately and strongly degenerate layers. The 
T s (Ti nt ) relation depends mainly on the conductivities in 
the sensitivity strip on the p — T plane (Gudmundsson et 
al. 1983) placed near the turning point, where n c ~ n r . 



2.3.1. Radiative opacities 

Radiative opacities of partially ionized iron in s trong mag- 
netic fields were studied by Rajagopal et al. ( 1997 ), but 
neither tables nor analytic fits required for NS modelling 
were published. Fortunately, we will see that the T s (Ti nt ) 
relation for a not too cold NS is not noticeably changed 
if we replace the radiative opacity of partially ionized 
iron by that of fully ionized iron. Therefore we base our 
opacity model on the assumption of full ionization. In 
this case, the two principal contributions into the opac- 
ity come from free- free transitions in electron-ion collisions 
and from Thomson scattering of photons by free electrons. 

At B = 0, the second process is described by the well- 
known Thomson cross section <tt in the non-relativistic 
limit T < T T (e.g., Berestetskii et al. |l992j ). The corre- 
sponding opacity is 



A" i 



ra c cr T 



8tt 

T 



(14) 
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The free-free absorption coefficient at B — was cal- 
culated by Karzas & Latter (1961) as a function of the 
electron velocity and photon frequency u). Hummer (1988) 
produced an accurate fit to its thermal average in the 
range 1(T 4 < £ < 10 15 and 1(T 3 < T Ry < 10 3 , where 



2U 2 



!Ry 



,Z 2 C 



kv,T 



Ta 



0.157 89 Z 2 ' 



TlLO 



(15) 



The ratio of the Rosseland mean, Kg, calculated using 
that fit, to the Thomson opacity can be written as 



Zotf 

3C7 



c- 



T^' (16) 



where p is measured in g cm , cvf — e / (he) ~ j^f, 



15 



(l-e-«)3A ff (T,0 



(17) 



Aff = (tt/^/3) gg is the Coulomb logarithm, and gg is the 
thermally averaged Gaunt factor. In the Born approxima- 
tion (e.g., Silant'ev & Yakovlev ^98^) , A ff = e^ 2 K (£/2), 
where Kq is a modified Bessel function; then C7 = 316.824. 
Since e^/ 2 K (£/2) w \/tt/£, at £ — > 00, we can obtain a 
reasonable approximation to gg by extending the accurate 
fit of Hummer ( |1988[ ) beyond its boundary £ = 10 1,5 in 
the power law form gg oc £ -1 / 2 . Then we can calculate 
the integral in Eq. @. The result is fitted (within 0.6%) 

by 



108.. 



C7 



77.6 T° y 834 



1 



0.502 T° y 355 



0.245 r° y 834 ' 



(18) 



If free-free transitions and scattering at B — are 
important simultaneously, the Rosseland mean of the sum 
of the two opacities can be presented as 

K r = (K s + K T )A(f,T), f = K s /(K s + K T ). (19) 

The non-additivity factor A(f, T) has been calculated by 
Silant'ev & Yakovlev (1980) in the Born approximation 
(in which case it does not depend on T) . Beyond the Born 
approximation, we obtain a function weakly depending on 
T, which can be fitted (within 1%) as 



A(f,T) = l 



1.097 + 0.777 T Ry 
1 + 0.536 T Ry 



/°- bl7 (W) 



0.77 



(20) 



Thus we can propose a simplified analytic treatment of 
the radiative opacity at B = as the opacity produced by 
free-free transitions and Thomson scattering alone. The 
effect of this simplified treatment is illustrated in Fig. |l|. 
Here, we have applied the code for solving the thermal 
structure equation, described in PCY, with the updated 
conductive opacities (Potekhin et al. 1999) and two dif- 
ferent sets of radiative opacities. The solid line on the 
upper panel is obtained using the OPAL radiative opac- 
ities (Iglesias & Rogers 1996 ), whereas the dashed line 
is obtained using Eqs. (|l4[)-(20). The resulting difference 
(lower panel) does not exceed several percent in the most 



6.5 



E-H 



E-H 

<l 



5.5 

0.1 


-0.1 




ff+Th. 



_L 



7 8 9 

ig r tot [k] 

Fig. 1. T s (T lnt ) relations at B = and g = 2.43 x 10 14 
cm s~ 2 , obtained using OPAL radiative opacities (solid 
line) and the free-free+Thomson opacities described in 
the text (dashed line). Lower window: fractional differ- 
ences between the accurate numerical values of T s (Ti nt ) 
and those obtained using (i) the free-free+Thomson opac- 
ity (dashed line) and (ii) the older analytic formula of 
PCY (dot-dashed line). 



important range of T; nt shown in the figure. The differ- 
ence between the curve obtained using Eqs. (|l4|)-(p0|) and 
the fit derived in PCY (dot-dashed line) is of similar mag- 
nitude. The latter difference arises not only from the fit 
error (PCY), but also from the improvement of the elec- 
tron conduction opacity (Potekhin et al. 1999] ) and, at 
T int ^ 10 9 K, from shifting p h to 4 x 10 11 gem" 3 (for, 
as noted in PCY, isothcrmality is not fully reached at 
p = 10 10 gcur 3 for high T). 

Interaction of photons with magnetized plasma de- 
pends on their polarization and propagation direction. 
Accordingly, radiative thermal conductivities and associ- 
ated opacities along and across the field become different. 
Rosseland mean opacities K r ^ and K r j_ were calculated by 
Silant'ev & Yakovlev ( 1980 ) for various values of T and 
B. We have fitted their numerical results by 



K l{l (p,T,B) 
K T (p,T,0) 
K lX (p,T,B) 



1 



Ai u + [A 2 uf 



1 + A3 u 2 

1 + (At u) 3 - 5 + (A 5 uf 
1 + A 6 u 2 



(21) 
(22) 



6 



A. Y. Potekhin and D. G. Yakovlev: Thermal structure and cooling of magnetic neutron stars 



Table 1. Coefficients a n , b n and c„ of Eq. (12 



n 


1 


2 


3 


4 


5 


6 


a„ 
b n 
C„ 


0.0949 
0.0610 
0.090 


0.1619 
0.1400 
0.0993 


0.2587 
0.1941 
0.0533 


0.3418 
0.0415 
2.15 


0.4760 
0.3115 
0.2377 


0.2533 
0.1547 
0.231 




A n 


= a n - 








(23) 



thermal conductivity is mainly determined by electron 
scattering off ions. At B = it can be written as 



where u = T B /(2T), K r (p,T,0) is the field-free opacity 
(|l9|), and the parameter / introduced in Eq. (|l9|) is cal- 
culated at B = 0. The best- fit parameters a n , b n and c n 
given in Table |l| ensure an average fit error of 5.5% with 
the maximum error of 11%. 

Asymptotic behaviour of Eqs. (21) and (E2) at u — > oo 
agrees with theoretical results of Silant'ev & Yakovlev 
j98C|) : K l{l (p,T,B) = 2K r± (p,T,B) = (ir/u) 2 K T (p, T, 0) 
at / = (Thomson scattering); K x \\(p, T, B) = 
K T ±(p,T,B) at / = 1 (free-free transitions); and K T cx 
u~ 2 at any /. 

At finite but large u the radiative opacities of fully ion- 
ized matter are strongly reduced. The reduction is ~ 10 
times stronger for scattering than for free- free transitions. 
Inserting this factor 10 into Eq. ([l6]) and taking into ac- 
count that, at B <; 10 11 G, the NS radiative sur face i s 
pushed to p s ~ B\2 gcm~ 3 (Ventura & Potekhin 2001 ), 
one can see that in deep, strongly magnetized photo- 
spheric layers Thomson scattering dominates the opacity 
only at T 6 > 10 p 2/T > 10 B% 7 . 

We note that the Rosseland mean opacity due to scat- 
tering by free ions, , can be obtained from that for 
electrons, Kt, by the simple scaling: 



K^(p,T,B)nZ 



3 / m e 



K T {p,T,—^B 
Am,, 



(24) 



On the contrary, there is no simple scaling for the opacity 
of photons due to electron-ion (electron bremsstrahlung) 
and ion-ion (ion bremsstrahlung) collisions. For instance, 
the processes are well known to be drastically different in 
the electric dipole approximation. 

From the asymptote of Kt at large u we see that the 
Thomson opacity of ions is ~ Z times larger than that of 
electrons, if the ions are strongly quantized into Landau 
levels, which occurs at T 6 < 0.07 (Z/A) B 12 . Since the 
main contribution into the radiative opacity at sufficiently 
low temperature comes from free-free processes (see 
above), we can conclude that the Thomson ion scatter- 
ing cannot contribute appreciably to the Rosseland mean 
opacity of the photosphere unless B\i ^> 10 3 (A/Z) 1A . 

2.3.2. Electron conductivities 

Electron thermal and electrical conduction is the most im- 
portant process that determines thermal structure and 
magnetic evolution of NSs. In dense, strongly coupled 
Coulomb plasmas, typical of NS envelopes, the electron 



3m* 



(25) 



yl + x 2 is the kinematic electron mass 



where m* = m, 
and t k is the effective relaxation time. 



Recently, Baiko et al. ( 1998 ) considerably improved 
the treatment of k by taking into account multiphonon 
absorption and emission processes in Coulomb crystals 
and incipient quasiordering in a Coulomb liquid of ions. 
Using these results, Potekhin et al. ( 1999 ) constructed 
an effective scattering potential which allowed them to 
calculate the non-magnetic conductivity in the relaxation 
time approximation in an analytic form; these analytic ex- 
pressions described accurately numerical results obtained 
beyond the framework of the relaxation time approxima- 
tion. The effective potential has been used then at ar- 
bitrary magnetic fields in order to derive practical ex- 
pressions for evaluation of electrical and thermal conduc- 
tivities of degenerate electrons in magnetized outer en- 
velopes of NSs (Potekhin 1999). Unlike previous treat- 



ments of the electrical conductivities perpendicular to the 



quantizing mag netic fie lds (K aminkcr & Yako vlev 1981 
Hernquist |l984j Schaaf |l988|) , Potekhin (|l99S|) went be 
yond the assumption that w g r K 3> 1 by introducing an 
interpolation from the case w g r K ^> 1 to w g r K -C 1. All 
tensor components of the kinetic coefficients in magnetic 
fields have been calculated and fitted by analytic formulae; 
the corresponding Fortran code is available electronically 
at http: / / www.ioffe.rssi.ru / astro / conduct / . This code is 
used here to calculate the temperature profile in the outer 
envelope of a NS. 

In the inner envelope, we need to know the thermal 
conductivity k only in young NSs, t ^ 100 yr, as long 
as the internal region (r < i?t>) remains non- isothermal 
(GYP). For this purpose, we will use the non-magnetic 
conductivities presented by GYP. They generalize the ex- 
pressions by Potekhin et al. (199£) in two respects. First, 
the atomic nuclei in the inner envelope cannot be con- 
sidered as pointlike: the size and shape of nuclear charge 
distribution significantly affect the conductivity. Second, 
the electron-phonon scattering in the inner crust at tem- 
peratures much below 10 s K changes its character: the so 
called umklapp processes cease to dominate and the nor- 
mal processes (with electron momentum transfer within 
one Brillouin zone) become more important. 

The internal temperature in young NSs is relatively 
high, T *Z 10 8 K. Our neglect of the effects of the magnetic 
fields in the internal region is strictly justified for the case 
of small magnetization parameter, uj g T K < 1. This param- 
eter is shown in Fig. |^ as a function of T for three values 
of density: near the top and the bottom of the inner enve- 
lope (solid lines) and in the middle of the outer envelope 
(dashed line). We see that the use of non-magnetic K is 
justified for the NSs with the magnetic fields B <J 10 13 G 
in the inner envelope. However it cannot be justified for 
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Fig. 2. Electron magnetization parameter at the top and 
bottom of the inner crust (solid lines) and in the outer 
envelope (dashed line). The curves are marked with the 
values of lgp [g cm -3 ]. 



higher B in the inner envelope or for B <; 10 10 G in the 
outer envelope. In such cases the conductivity across the 
field is suppressed by the large factor [1 + (ui g T K ) 2 ] and the 
Hall thermal conductivity may become important. We ne- 
glect these effects in the inner envelope. Therefore, our 
simulations (Sect. 4) of cooling of young NSs [t & 100 
yr) are justified for B 10 13 G in the inner envelopes. 
The results for older NSs, for which the inner envelope is 
isothermal, are valid at larger B as well. 

3. Thermal structure 

3.1. Temperature profiles 

We have solved Eq. (Q) using the numerical algorithm de- 
scribed in PCY with the EOS and opacities presented in 
the previous section. We have varied the magnetic field in- 
clination angle 9 from to 90°, the effective temperature 
T s from 2 x 10 5 K to 10 7 K, and the field strength B from 
to 10 16 G. As clear from the discussion in Sect. 0, our 
model EOS and opacities may be crude at B > 1CF 4 G. 
Further improvements of physics input are required for 
obtaining more reliable results at such B. 

Figures || and |] present the calculated T(p) profiles 
for a NS of mass M = 1.4M Q and radius R = 10 km 



(g = 2 A3 x 10 1 



and r g /R — 0.413). The curves in 



Fig. U are calculated at fixed T s = 10 6 K. The solid curves 
show the results of accurate calculations for the magnetic 
field normal and tangential to the surface (cases of parallel 
and transverse conduction, Sect. 2.1.2). The dashed lines 



are obtained for the parallel conduction using the classi- 
cal thermal conductivity (neglecting the Landau quanti- 
zation). This approximation becomes inaccurate with in- 
creasing B. The dot-and-dash lines are calculated for the 
tangential magnetic field using the transverse conductiv- 
ity in the limit o f oj z t k ;g> 1, employe d in previous papers 



(e.g., Hcrnquist |l984j |1985|; S chaaf |l988l |l990a| , |1990b 
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Fig. 3. Temperature profiles through an iron envelope of 
a NS with M = 1.4M and R — 10 km for the magnetic 
field B = 10 11 G (left panel) and 10 13 G (right panel) 
at the fixed effective surface temperature T s = 10 6 K. 
Lower solid curves correspond to 9 = and upper ones 
to 9 = 90°. Dashed lines are calculated with the classical 
thermal conductivity for 9 = 0, while dot-dashed ones 
are calculated for 9 = 90° assuming oj g r K ^> 1 which is 
traditional for the transverse conduction. 



line) is given by the interpolation (Potekhin 1999 ) men- 
tioned in Sect. |2.3.2 . 

Figure [| shows temperature profiles at four field 
strengths B and five inclinations 9 for T mt = 10 8 and 
10 9 K (we do not present the curves for T; nt = 10 8 at 
B = 10 15 G because the temperature T(p) is too small 
near the surface in the transverse conduction case). The 
curves show pronounced ^-dependence at high densities 
starting from the turning point (Sect. pTS ) . 

The endpoints of the curves in Figs7"| and ^ lie at the 
radiative surface p = p s , where T = T s , and near the neu- 
tron drip point. Their behaviour qualitatively agrees with 
the results o f the approximate analytic study of Ventura 
& Potekhin ( 2001 ). For instance, we confirm the approxi- 
mate linear dependence p s cx B and the shift of the turning 
point to higher densities with increasing B. The density 
Pb marked on the graphs is defined by Eq. (|]). It sep- 
arates the strongly- and weakly-quantizing field regimes. 
We see that the i ncreas e of T at p > pb (neglected by 
Heyl & Hcrnquist 1998a ) can be significant. 

The higher the temperature, the wider is the den- 
sity region, where k\\ and kj_ are of similar magnitude. 
Therefore the dependence of the profiles on the inclina- 
tion 9 and magnetic field strength B is less pronounced at 
higher T mt showing convergence to the B — case. 



3.2. Relation between internal and effective 
temperatures 

A simple analytic fit t o T s (T m t) 



was 



Heyl & Hernquist 1998a , 1998b ) . This approximation is in- 
accurate at lower B; a more accurate approximation (solid 



found by 

Gudmundsson et al. (1982) for iron envelopes at T s > 
3 x 10 5 K. Later, using an updated physics input, PCY 
constructed a more general fit for NSs with (accreted) en- 
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Fig. 4. Temperature profiles through an iron envelope of a NS with M = 1.4 Mq and R — 10 km at (from left to right) 
B = 10 u , 10 12 , 10 14 , and 10 15 G. The internal temperature is fixed to T int = 10 8 K (solid lines) or 10 9 K (dot-dashed 
lines). The lines of each group correspond to cos# = 1 (the lowest line), 0.7, 0.4, 0.1, and (the highest line). 



velopes containing layers of different chemical elements. 
In all cases, the scaling law T s cx g 1 ^ 4 has been confirmed. 
This approximate law, first noted by Gudmundsson et al. 
( |1983| ) , follows from the thermal structure equation , if 
one takes into account that the largest contribution to the 
integral of this equation comes from the sensitivity strip, 
being almost independent o f the outer boundary condi- 
tion (cf. Ventura & Potekhin |oo|) . The scaling T s cx g 1 / 4 
holds also in our case of magnetized envelopes. Therefore 
the ratio X of the surface temperature T S (B, 9, <?, Tj nt ) at 
a given magnetic field to the value Ts°\g, T- mt ) in the ab- 
sence of the field is practically independent of g: 



T S {B, 9, g, T int ) « Ti%, T int ) X(B, 6, T int ). 



(26) 



Since our outer envelope consists of iron, 

T (0) 

is given by 

Eq. (A7) of PCY: 



rw^io 6 ^ 4 [(7C) 2 - 25 + (C/3) x 



.251 V4 



K. 



(27) 



where ( = T int , 9 - 0.001 g{{ 4 y/TF^, T int , 9 = 
T int /(10 9 K) and g u = g/10 14 cm s" 2 . 

For the parallel and transverse conduction cases (at 
the magnetic pole and equator) we have, respectively, 



Xu(B,T iz 



1 + 0.0492 B° 12 292 /T^j°, 



(28) 



v (R T . \ ~ + 0-1076 Bi 2 (0.03 + r int , 9 )-°-^ 
±{ ' mtJ ~ [1 + 0.819 5 12 /(0.03 + T int , 9 )]°- 6463 ' l ' 

These fits have been checked for B < 10 16 G and 10 7 K < 
Tmt < 10 95 K with an additional constraint T s > 2 x 
10 5 K. Maximum residuals (11% and 5%, respectively) 
occur at the largest B. 

The effect of the magnetic field on the effective tem- 
perature and the quality of the present fit are illustrated in 
Fig. ||. The magnetic field normal to the surface enhances 
the photon luminosity, since heat is transported along the 
field lines (||). Accordingly, T s grows up with increasing B. 
In contrast, the tangential magnetic field significantly de- 
creases T s since the heat propagates across the field lines 



(_L). Our results arc in qualitative agreement with earlier 
results of Schaaf ( |1990aD and Heyl & Hernquist ( |1998b[ ), 
also shown in Fig. |5|, although there are quantitative dif- 
ferences. The origin of the large difference from Schaaf 



(1990a) is not clear, whereas the difference from the re- 
sults of Heyl & Hernquist ( 1998b| ) can be attributed to our 
more accurate treatment of electron thermal conductivity, 
particularly in the nondegenerate region of the envelope. 
Note that the results of Heyl & Hernquist ( 1998b ) for par- 
allel conduction do not converge to the B = results in 
the limit of high Tint - In addition to the above authors, the 
case of par allel c onduction was considered thoroughly by 
Van Riper (198S) who however did not present any tables 
or fit expressions for practical applications. We have com- 
pared our results with those presented in his Fig. 29. In a 
hot NS (Tint ^ 3x 10 9 K) his results are in good agreement 
with ours and satisfy the criterion of convergence to the 
B = case. For B = the agreement is also quite good 
at any T int £ 10 7 K (T s £ 3 x 10 5 K) but for B £ 10 13 
G and Tint C3x 10 9 K the magnetized envelopes of Van 
Riper appear to be much more heat transparent than ours 
(in contrast, the envelopes of Heyl & Hernquist 1998b| are 
less transparent than ours at the same conditions). The 
nature of this difference is also unclear. 



3.3. Variation of temperature over stellar surface 

The dependence of T s on the angle 8 between the magnetic 
field and the normal to the surface is most easily described 
by the model of Greenstein & Hartke (1983), which im- 
plies a superposition of "longitudinal" and "transverse" 
heat fluxes: T*(6) = T 4 (0) cos 2 6> + !T 4 (90° ) sin 2 9. This ap- 
proximation has been used, e.g., by Page (|1995 ); Shibanov 
& Yakovlev fll996j ); Heyl & Hernquist ( |1998aj ). Our nu- 
merical calculations confirm that it reproduces accurately 
(within w 30%) the T s (6) dependence. However, a replace- 
ment of the power index 4 with 4.5 yields better accuracy 
(within 10%): 

X{B,6,T- mt ) = [X* /2 (B,T int )cos 2 
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Fig. 5. Dependence of the surface tem- 
perature T s on the internal tempera- 
ture Tint at the magnetic fields B = 0, 
10 13 , and 10 14 G normal (left panel) 
or tangential (right panel) to the sur- 
face. Numerical data (heavy dots) are 
compared with the present fit (pri|)-(|29|) 
(solid line s) and with the earlier fits 
of Schaaf (1990a) (dot-dashed line, for 



B = 1 0" G) and Heyl & Hernquist 



(1998b) (dashed lines). 



Xl /2 (B, T int ) sin 2 



1 2 /9 



(30) 



Still better accuracy can be reached by varying the power 
index between 4 and 5 as a function of B and T; nt . 

Figure || shows the distribution of the effective tem- 
perature from the magnetic pole to the equator. Heavy 
dots show the results of numerical integration of T(p) pro- 
files. Solid lines are obtained using our fit, Eqs. (]2q)-(p0|), 
while short-dashed lines correspond to the approximation 
of Greenstein-Hartke ( |1983 ). For comparison, dot-dashed 
lines on the left panel represent the fitting formulae of 
Schaaf ( pJ90a| [l990b| ). 

Horizontal long-dashed lines in Fig. || show T s which 
would have been observed at the same Tint if the mag- 
netic field were absent. In qualitative agreement with the 
results of earlier works (Sect. Q), the thermal flux emer- 
gent from the NS interior is suppressed by the magnetic 
fields in the equatorial region and enhanced near the pole. 
For a rotating NS, su ch dis tribution leads to a pulsating 
light curve (e.g., Page 1995 ), which may be observed, for 
instance, with Chandra in soft X-rays (where the spectral 
flux of the thermal NS radiation has maximum) or with 
the HST in UV. 



3.4. Total photon luminosities 

In order to calculate the cooling curves, one has to know 
the total NS photon luminosity L at any given Tj nt . In 
the non-magnetic case, the radiation flux emergent from 
the surface is given by F = L/(4irR 2 ) = <tT s 4 , where T s = 
Ts°\g, Ti n t) is known from Eq. (|27j). In the magnetic case, 
F should be averaged over the surface. Using the fitting 
formulae (p8|)-(30), we have performed such averaging for 



the dipole configuration of the field [Eq. (|i"fj|)1 at B p < 
10 16 G, T int = (10 7 -10 9 5 ) K, and r g /R < 0.7. The ratio 
of magnetic to non-magnetic fluxes can be fitted as 



F(B) _ 1 + ai/3 2 + a 2 f3 3 + 0.007 a 3 (3 4 
where 



1 + 



(31) 



= 0.074 y/B^Tr°f, 



5059 T 



.3/4 



ai 



«2 



«3 



int. 9 



(1 + 20.4T; 



.1/2 



int,9 



138 T 



,3/2 



int. 9 



•11022* )V* 



1484 T 



,3/4 



(1 + 90^4 
5530(1 



125 



)l/2 

0Ar g /R)-^T^ 9 



(1 + 8.16 2T 1 n f 9 + 107.8 T. 



n3/2 



Mnt,9 



-560T 2 tj9 )V2 



and B\2 relates to the magnetic pole. The maximum error 
of this fit is 6.1% (i.e., 1.5% for T ). Note that Eqs. @- 
( |30| ) used for calculating F(B) are actually less accurate 
by themselves, which lowers the real accuracy of our an- 
alytic fits. Nevertheless the latter accuracy seems to be 
sufficient for cooling simulations. 

Note one important feature: the effect of the magnetic 
field on the F(B) /F(0) ratio becomes weaker with growing 
Tjnt ■ It is explained by the arguments presented in Sect. 
|3.l[ Accordingly, the luminosity of a hot NS cannot be 
strongly affected even by very high magnetic fields. 

4. Cooling 

Let us discuss briefly the effects of magnetized envelopes 
on NS cooling. For illustration, we take the models of NSs 
with two masses, M — 1.3 and 1.5 Mq. For the chosen 
EOS in the stellar core (Sect. |), the 1.3 Mq NS has the 
radius R — 11.86 km. Its central density p c — 1.070 x 10 15 
g cm -3 is lower than the density pdu = 1-298 x 10 15 g 
cm~ 3 at which the powerful direct Urea process of neu- 
trino emission is switched on. Thus the main neutrino 
emission mechanisms in the NS core are a modified Urea 
process and neutrino bremsstrahlung in nucleon-nucleon 
collisions. Accordingly, the 1.3 Mq model gives us the ex- 
ample of slow cooling. For M — 1.5 Mq we have R = 11.38 
km and p c = 1.420 x 10 15 g cm -3 > pdu- The direct Urea 
process is open in the central kernel of mass 0.065 Mq and 
radius 2.84 km. This is an example of fast cooling. Further 
details of these NS models may be found in GYP. 

We have calculated the cooling curves of NSs with 
magnetized envelopes, assuming the dipole magnetic field, 
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Fig. 6. Dependence of the surface tem- 
perature T s on the inclination an- 
gle 6 at constant B and Tint- Left 
panel: B = 10 12 G, lg T int [K] = 8.0, 
8.5, 8.8; right panel: B = 10 14 G, 
lg Tint [K] = 8.7. Dots: numerical re- 
sults; solid lines: fit; dashed lines: in- 
terpolation of Greenstein-Hartke be- 
tween the accurate values at 9 — 
and 90°; long dashes: the case of B — 
0; dot-dashed lines on the left panel: 
combined fitting formulae of Schaaf 
(1990a,1990b). 



Eq. (|io|). For comparison with earlier papers (Sect. [T]), we 
also consider a familiar model of the field which is radial 
and constant in magnitude throughout the blanketing en- 
velope. In addition, we have switched on and off superflu- 
idity of neutrons and protons in stellar interiors to demon- 
strate the combined effects of magnetized envelopes and 
superfluid interiors on NS cooling. 



4.1. Cooling of non-superfluid neutron stars 
4.1.1. General features 

Figure ^| shows the cooling curves for non-superfluid 
1.3 Mq and 1.5 Mq NSs with dipole magnetic fields. The 
curves are marked with the magnetic field strengths at the 
magnetic poles, B pi and display the effective temperature 
r e °°, Eq. (g). During the first 50 yr of their lives, 1.3 Mq 
and 1.5 Mq NSs have nearly the same surface temper- 
atures since the stellar interiors are non-isothermal and 
thermal evolution of the stellar crust is decoupled from 
the evolution of the core (e.g., GYP). Later, thermal equi- 
librium is established in the interior, and the surface tem- 
perature of the 1.5 Mq NS becomes much lower due to 
very powerful direct Urea processes in the stellar kernel. 

The results for B < 10 14 G are in good qualita- 
tive agreement with those obtained by Page (1995) and 
Shibanov & Yakovlev (1996). The thermal state of the 
stellar interior is almost independent of the magnetic 
field in the NS envelope at the neutrino cooling stage 
(t ^ 10 4 — 10 5 yr), but is affected by the magnetic field 
later, at the photon cooling stage. On the contrary, the 
surface temperature is always affected by the magnetic 
field. The dipole field B ^ 10 13 G makes the blanketing 
envelope overall less heat-transparent. This lowers T c at 
the neutrino cooling stage and slows down cooling at the 
photon cooling stage. The dipole field B ^> 10 13 G makes 
the blanketing envelope overall more transparent, increas- 
ing T e at the neutrino cooling stage and accelerating the 
cooling at the photon stage. 




t (yr) 

Fig. 7. Cooling curves of non-superfluid 1.3 Mq and 

1.5 Mq NSs with dipole magnetic fields of different 
strengths. 



Contrary to the case of the dipole magnetic field, any 
radial field (e.g., Van Riper |1991| ) would always lower the 
thermal insulation of the blanketing envelope, increasing 
T s at the neutrino cooling stage and accelerating cooling 
at the photon cooling stage. The radial field would affect 
the cooling noticeably more than the dipole field. 

Even very strong magnetic fields do not change sig- 
nificantly thermal insulation of a hot blanketing envelope 
(Sect. 3J5). Accordingly, the effects of magnetic fields in a 
young and hot NSs are not too strong. The strongest ef- 
fects take place in cold NSs, at the photon cooling stage. 
Unfortunately, our knowledge of insulating properties of 
the blanketing envelopes is the poorest for these NSs (Sect. 
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g). In addition, the magnetic field evolution and reheating 
processes may become important for these stars, which 
we ignore for simplicity. Notice that the fields B <; 10 14 G 
appreciably accelerate the cooling and lower the duration 
of the neutrino cooling stage. 



4.1.2. Effects of magnetic field and light bending 

Let us consider some effects of the magnetic fields on NS 
cooling in more detail. Figure |^ demonstrates the effect 
of dipole or radial magnetic fields on the flux of electro- 
magnetic radiation F(B) detected from a non-superfluid 
1.3 Mq NS of an age of 25 000 yr. This is the age of the 



Vela pulsar (Lyne et al. 1996). This is also a typical age of 
anomalous X-ray pulsars which are likely to be magnetars 
(e.g., Mereghetti 2001). Any radial magnetic field is seen 
to increase the flux. For the field strength B = 4 x 10 12 G, 
which can be appropriate for the Vela pulsar, the increase 
is by a factor of 1.5, while for ultrahigh B = 10 15 and 10 16 
G the factor is about 4.5 and 9, respectively. 

For the dipole field, however, the detected flux depends 
on observation direction. The dipole field B p = 4 x 10 12 
G reduces the flux averaged over the entire NS surface 
(solid line) by a factor of 1.3, while the fields B p = 10 15 
and 10 16 G amplify it by factors of about 2 and 4, respec- 
tively. The dashed and dotted lines in Fig. || refer to the 
fluxes detected in the direction to the magnetic pole and 
equator, respectively (assuming the magnetic axis coin- 
cides with the rotational one, for simplicity). These lines 
show the largest difference of the fluxes detected under dif- 
ferent angles. They are obtained taking proper account of 
the general relativity effect of bending of light rays prop- 
agating from the NS surface to a distant observer, in the 



same manner as was done by Page (1995). For comparison, 
we present also the fluxes calculated neglecting the grav- 
itational ray-bending effect (as if space-time outside the 
NS were flat). In the absence of light bending effects, the 
difference of fluxes observed under different angles is quite 
noticeable. For instance, at B p = 4 x 10 12 G the largest 
difference would be about 65%. The light bending reduces 
the largest difference to 14% for B p = 4 x 10 12 G, mak- 
ing it practically negligible. This is explained (e.g., Page 
1995| ) by the fact that gravitational bending of light rays 



increases the fraction of the NS surface visible by the ob- 
server. The observer detects the flux from the larger por- 
tion of the surface, that is closer to the flux averaged over 
observation directions. The light bending effect is strong 
for our 1.3 Mg and 1.5 Mq NS models. Thus we will ne- 
glect weak dependence of fluxes on the detection direction 
and use the average fluxes and associated effective temper- 
atures T c in our analysis. Notice that here we mean total 
(spectral integrated) fluxes but not the fluxes observed 
in selected spectral bands taking into account interstellar 
absorption. 

Figure || shows the effective temperatures of young 
(1000 yr) non-superfluid 1.3 Mq and 1.5M NSs with 
dipole and radial magnetic fields of different strengths. 



1.0 



0.8 



o 0.6 

Si. 

§0.4 
1=0 

^ 0.2 



0.0 = 



-0.2 



1.3 Mq t=25000yr 



surf, aver 



Icnsing outside NS 



no lensm 




. ' " pole obs. 



equator obs. 



10 11 12 13 14 15 

lgB[G] 



16 



Fig. 8. Effect of dipole or radial surface magnetic fields on 
the flux of electromagnetic radiation F(B) detected from 
a non-superfluid 1.3 Mq NS of age t = 25000 yr. Solid 
line shows the flux averaged over the entire NS surface 
for the dipole field. Dot-and-dashed lines present the flux 
observed in the direction towards the magnetic pole (as- 
sumed to be aligned with the rotational axis) or towards 
the magnetic equator. The dotted lines show the same 
fluxes calculated neglecting gravitational bending of light 
rays outside the NS. 



The effects of magnetic fields are similar to those in Fig. 
||. Since the 1.5 Mq NS undergoes fast cooling, its in- 
ternal temperature (~ 10 7 K) is much lower than in 
the 1.3 Mq star (~ 3 x 10 s K). Accordingly, the mag- 
netic fields affect the surface temperature of the 1.5 Mq 
NS more significa ntly. Pr eviously, it was suggested (e.g., 
Heyl & Hernquist 1997a ) that the models of young and 
hot, slowly cooling NSs with ultramagnetized envelopes 
could serve as models of anomalous X-ray pulsars and soft 
gamma repeaters. Although those previous cooling models 
were rather simplified, they are in qualitative agreement 
with our more elaborate models. The main outcome of 
these studies is that even ultrahigh magnetic fields cannot 
change appreciably the average surface temperatures of 
young NSs with iron envelopes (as seen clearly from Figs. 
U and ||). We stress that the local surface temperatures 
in narrow strips near the magnetic equators appear to be 
much lower than near the magnetic poles (Sect. ||), but 
they do not contribute noticeably to the NS luminosity 
integrated over the surface. 
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Fig. 9. Effect of dipole or radial surface magnetic fields 
on the averaged effective surface temperature r o °° of non- 
superfiuid 1.3 Mq and 1.5 M NSs of age t = 1000 yr. 



4.2. Effects of superfluidity and application to the Vela 
pulsar 

Our analysis shows that the magnetic fields of ordinary 
pulsars (B ^5 10 13 G) do not strongly affect the sur- 
face temperature, at least at the neutrino cooling stage. 
Nevertheless, such fields may influence the theoretical in- 
terpretation of observations. For illustration, let us analyse 
observations of the thermal radiation from the Vela pul- 
sar. We adopt the value of the effective surface tempera- 
ture r o °° = (6.8 ±0.4) x 10 5 K inferred most recently from 
observations with the Chandra observatory by Pavlov et 



al. (2001) (at the la level) using the hydrogen atmosphere 
models. This value of T e °° is in a good agreement with the 
value T„°° = (7.85 ±0.25) x 10 5 K obtained earlier by Page 



et al. ( 1996 ) from ROSAT observations. Thus we assume 
the presence of a very thin hydrogen atmosphere (of mass 
^5 10~ 13 Mq) which determines the spectrum of thermal 
radiation from the NS surface but does not affect ther- 
mal insulation of a blanketing envelope composed mostly 
of iron. Let us use the slowly cooling 1.3 Mq NS model 
for the interpretation of observations. Assuming B = 
and no superfluidity in the NS interior, we have, for the 
age of 25 000 yrj]!; 00 » 1.0 X 10 6 K, which is higher than 
the value obtained from observations (Fig. [l(]) . The dipole 
magnetic field B p = 4x 10 12 G lowers T c °° by 8% (Fig. g), 
not sufficient to attain the required value. 



1 We have verified that T e °° is almost independent of the NS 
mass, for the given EOS, as long as M is lower than the critical 
value 1.44 Mq at which the direct Urea process becomes open 
in the stellar interior. 
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Fig. 10. Cooling of 1.3 Mq NS with B = 0, with the dipole 
magnetic field B p = 4 x 10 12 G or with the radial field of 
the same strength. The upper curves are calculated as- 
suming no superfluidity in the NS interior while the lower 
ones are calculated for certain critical temperatures of the 
neutron and proton superfluids in the NS cores (the val- 
ues of T cn and T cp are indicated near the curves in units of 
10 8 K). Error bars show the possible interval of the surface 
temperatures of the Vela pulsar (see text). 



It is well known that neutron or proton superfluidity in 
the NS core affects the cooling. For illustration, we adopt 
the same simplified model of nucleon superfluidity which 
ha s been used in a number of previous works (Yakovlev et 
al. 1999 and references therein). In this model, the neutron 
pairing is assumed to be in the triplet state, while the pro- 
ton pairing occurs in the singlet state of a nucleon-nucleon 
pair. The critical temperatures T cn and T cp of the neutron 



and proton superfluids are assumed to be constant and 
treated as free parameters to be adjusted to observational 
data. Figure [l0| demonstrates that taking the neutron and 
proton superfluidities with certain values of T cn and T cp 
we can lower T c of our Vela model to the observed val- 
ues. This can be done for all three models of the envelope, 
with B = as well as with the dipole and radial fields. 
In the given examples, the lowering of T e for the dipole 
or field-free cases is produced by neutrino emission due to 
Cooper pairing of neutrons, while for the radial field the 
lowering is produced by neutrino emission due to Cooper 
pairing of neutrons and protons. In these three cases we 
need different critical temperatures T cn and T cp . In each 
case the choice of T cn and T cp is not unique. 

Generally if we fix the blanketing envelope model, 
we can determine the domains of T cn and T cp values in 
the T cn -T cp plane which force the NS to have T e °° within 
the given errorbar by the specified age t. These domains 
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Fig. 11. Domains (hatched) of T^ 
1.3 Mq NS with the superfluid core has the surface tem- 
perature r o °° = (6.8 ± 0.4) x 10 5 K by the age of the Vela 
pulsar. Solid lines enclose the domains for the NS with 
the dipole magnetic field (B p = 4 x 10 12 G); dashed lines 
are for the non-magnetized NS and dots are for the radial 
field {B = B p ). 



(hatched regions) are shown in Fig. |Tl| for the Vela pulsar 
with a non-magnetic envelope, and also for the envelope 
with dipole and radial magnetic fields. 

The domain for the non-magnetic envelope has a com- 
plicated X-shape. The dipole magnetic field facilitates a 
lowering of T c to the required errorbar. The domain of ac- 
ceptable values of T cn and T cp becomes narrower overall 
and splits into two parts (Fig. [ll]). In principle, in these 
two cases we do not need proton superfluidity to explain 
the observations (Fig. |l^): the neutrino emission due to 
Cooper pairing of neutrons is sufficiently strong by itself 
to lower T e . On the other hand, were the magnetic field 
radial, it would raise T e and complicate the adjusting of T e 
to observations. Neutrino emission due to Cooper pairing 
of neutrons is insufficient to lower T , and the emission 
due to Cooper pairing of protons is needed. The domain 
acquires a qualitatively different shape: the proton super- 
fluidity cannot be too weak (lgTcp ^ 8.3). 

The domains of T cn and T cp obtained for the Vela pul- 
sar can be combined further with analogous domains de- 
termined for other cooling NSs in attempt to impose re- 
strictions on the values of T cn and T cp in the stellar cores. 
A preliminary a nalys is of such a kind has been done by 
Yakovlev et al. (|l999|) . However, further work remains to 
be done to constrain strongly the critical temperatures: 
one should take into account new observational data on 
thermal radiation from isolated NSs and new theoretical 
models of cooling NSs with superfluid cores; particularly, 



one has to incorporate density dependence of T cn and T cp 
over the stellar core. Such work is beyond the scope of 
the present paper, but our illustrative examples show that 
proper analysis of the NS superfluidity may require con- 
sideration of the effects of the magnetic field on the NS 
blanketing envelope. 



5. Conclusions 

We have analysed the thermal structure of neutron star 
envelopes with typical pulsar magnetic fields 10 11 -10 13 G 
and with ultrahigh magnetic fields up to 10 16 G. We have 
used (Sect. ||) modern data on equation of state and ther- 
mal conductivities of magnetized neutron star envelopes. 
In particular, we have proposed an analytic model of the 
radiative thermal conductivity limited by the Thomson 
scattering and free-free absorption of photons in a magne- 
tized plasma. We have used the values of thermal conduc- 
tivity of degenerate electrons, updated recently by proper 
treatment of dynamical ion-ion correlations which affect 
electron-ion scattering. We have calculated the temper- 
ature profiles (Sect. ||) in the neutron star envelope for 
any inclination of the magnetic field to the surface and 
obtained a fit expression which relates the internal tem- 
perature of the neutron star to the local effective surface 
temperature. We have also calculated and fitted the rela- 
tion between the internal temperature and total surface 
luminosity (or effective temperature) for the dipole mag- 
netic fields B <, 10 16 G. 

Furthermore, we have performed (Sect. ^) simulations 
of cooling of neutron stars with dipole or radial magnetic 
fields in the envelopes. In agreement with the previous 
studies of Page ( 1995 ) and Shibanov & Yakovlev ( 1996 ), 
we have found that the effects of the two magnetic field 
configurations on neutron star cooling are qualitatively 
different. We have briefly discussed cooling of young and 
middle-aged neutron stars with ultramagnetized envelopes 
(magnetars as the models of soft gamma repeaters and 
anomalous X-ray pulsars) and the effect of surface mag- 
netic fields on constraining critical temperatures of the 
neutron and proton superfluids in the cores of ordinary 
pulsars, considering the Vela pulsar as an example. 

We stress that our results are less reliable for cold neu- 
tron stars (T s ^ 3 x 10 5 K) with very strong magnetic 
fields B <; 10 13 G because of lack of knowledge of ion- 
ization state and thermal conductivity of the outermost 
parts of the cold blanketing envelopes. Further work is 
required to fill these gaps. In this paper, we have consid- 
ered the blanketing envelope composed of iron and have 
not analysed the envelopes containing light elements. The 
presence of light elements generally reduces thermal in- 
sulation of the envelope. This effect has been studied in 
detail for non-magnetic envelopes (PCY) and is expected 
to be imp ortant for magnetic envelopes as well (Heyl & 
Hcrnquist 1997a ) . We are planning to address this matter 
in a future work. 
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